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We study a strongly correlated fermionic model with attractive interactions in the presence of disorder in two 
spatial dimensions. Our model has been designed so that it can be solved using the recently discovered meron- 
cluster approach. Although the model is unconventional it has the same symmetries as the Hubbard model. 
Since the naive algorithm is inefficient, we develop a new algorithm by combining the meron-cluster technique 
with the directed-loop update. This combination allows us to compute the pair susceptibility and the winding 
number susceptibility accurately. We find that the s-wave superconductivity, present in the clean model, does 
not disappear until the disorder reaches a temperature dependent critical strength. The critical behavior as a 
function of disorder close to the phase transition belongs to the Berezinky-Kosterlitz-Thouless universality class 
as expected. The fermionic degrees of freedom, although present, do not appear to play an important role near 
the phase transition. 

PACS numbers: 74.78.-w, 71.10.Fd, 02.70.Ss 
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I. INTRODUCTION 

A variety of systems show quantum coherence over large 
distances at low temperatures. Superfluidity and supercon- 
ductivity are two striking physical phenomena showing such 
behavior, which have been extensively studied over the years. 
However, when correlations between the microscopic degrees 
of freedom become strong it is difficult to study these phenom- 
ena theoretically from first principles. The calculations must 
take into account strong fluctuations over many length scales 
which is only possible numerically. When the microscopic 
degrees of freedom involve bosonic variables one can usu- 
ally devise efficient quantum Monte Carlo methods to solve 
the problem. 1 On the other hand, it is still difficult to study a 
variety of models from first principles when the microscopic 
theory is fermionic. For example, the critical temperature be- 
low which superconductivity is seen in the well-known at- 
tractive Hubbard model was only determined recently, 2 using 
the determinantal Monte Carlo method 3 - 4 and on lattices only 
as large as 18 x 18. The main approaches to dealing with 
fermionic systems can be viewed as arguments that univer- 
sality allows one to replace the microscopic theory with an 
effective low energy theory. The resulting effective theory is 
usually either a Fermi-liquid theory, a BCS-type mean-field 
theory, or some bosonic theory. 5 A key element in further- 
ing microscopic understanding, then, is to validate the univer- 
sality arguments and determine the low-energy effective the- 
ory; in practice, this has proved very difficult for systems with 
strong correlations. 

Real systems usually contain impurities. Thus, in addition 
to understanding superconductivity in clean systems, the ef- 
fects of impurities in the form of disorder need to be incorpo- 
rated in the studies. In certain systems like two-dimensional 
superconducting films and Josephson-j unction arrays, it has 
been discovered that superconductivity can be destroyed by 
tuning parameters such as the film thickness. 6 ' 7 Since these 
tuning parameters change the effective strength of the dis- 
order, it is believed that the superconductor-to-insulator (SI) 
phase transitions in these systems can be understood as be- 
ing driven by disorder. Among the models used to explain 



the experiments, the attractive Hubbard model with disordered 
chemical potential is one well-known starting point. 8,9 

The relevance of disorder for superconductivity was first 
addressed by Anderson, 10 where he argued that superconduc- 
tivity is insensitive to perturbations that do not destroy time 
reversal invariance. Using a BCS type trial wave function Ma 
and Lee 11 showed that superconductivity can persist even be- 
low the mobility edge. Clearly, these studies suggest that an SI 
transition is an effect of strong disorder which makes it a diffi- 
cult subject for analytic work. Fisher et alM. have argued that 
the effective theory describing the transition is bosonic, and 
then developed a deeper understanding of the purely bosonic 
superfluid-insulator transition using renormalization group ar- 
guments along with scaling. A variety of quantum Monte 
Carlo work has been done over the years on these purely 
bosonic microscopic theoriesii^i 4 *^ 6 *!^^ If fermions do not 
play an important role near the transition, it is likely that these 
studies will also be useful in understanding the universality 
of the fermionic SI transition. Recently determinantal quan- 
tum Monte Carlo studies of the attractive (fermionic) Hubbard 
model with disorder have been performed, 9 which show that 
it is indeed possible to drive an SI transition by increasing 
the disorder and, as expected, the critical disorder is large. 
However, the system sizes explored were quite small, 8x8. 
Other studies of disorder effects also involved only very small 
systemsjA*^ 

Motivated by the physics of the SI transition, in this arti- 
cle we study the effects of disorder in a strongly correlated 
fermionic model. Our model is unconventional and has been 
built so that it can be studied using the recently discovered 
meron cluster algorithms. 21,22 These novel algorithms are so 
efficient that lattices as large as L = 128 were studied recently, 
and it was shown with great precision both that the low tem- 
perature phase of the clean model is indeed superconducting 
and that the finite temperature phase transition belongs to the 
Berezinski-Kosterlitz-Thouless (BKT) universality class^ 3 ^ 
Here we use the same model to explore the effects of disor- 
der on superconductivity and, in particular, focus on the role 
of fermions. Unfortunately, the naive extension of the ear- 
lier algorithm becomes inefficient in the presence of disorder; 
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hence, we also develop a new algorithm by combining the 
meron-cluster formulation with the directed loop algorithm. 25 
This new algorithm allows us to measure the relevant observ- 
ables very accurately. 

Our paper is organized as follows: In Section II, we intro- 
duce our model and define the observables that we use later. 
In Section EH, we rewrite the model in a cluster representa- 
tion. Section IV explains the new directed-loop algorithm we 
have developed. Section V contains our results, and Section 
VI contains our conclusions and directions for the future. 



II. THE MODEL 

The model we consider in this article was motivated by 
the ability to solve the fermion sign problem using the the 
meron-cluster algorithm. 22 The Hamiltonian of the model can 
be written as 

(*j> » 

(2) 

where H-j consists of all the nearest neighbor interactions 
between sites i and j on an L x L square lattice and 
includes interactions on the site i. The term is given by 

= (^+J r 3-l)(n, T -i)(n ii -i) 

where U represents the Hubbard interaction between spin-up 
and spin-down electrons and /Ltj is the local chemical potential, 
through which disorder is introduced in the model. The term 

(2) 

H^j is unconventional and is given by 

H\? = i(cJ tT c i£r +ct CT c iCT ){(l + J 3 ) (4-47^+3) 
-(1- J 3 )K -2)} 
+ {Si-S^+Ji-Jj-Cl- J 3 )JfJf} (3) 

Here c ia and Ci a are the usual creation and destruction opera- 
tors of spin a at site i,n= c'c, and n,j = n i" + n j<?- 
is the spin operator on site i defined by 

Si = ~^24 s a ss ,c is < (4) 

s,s' 

and is the pseudo-spin operator defined by 

= {-iy* +3 *c l[Clh (5) 

J l = ^(^T +mi - !)• 



J + and J~ are related to pair creation and annihilation opera- 
tors. In our notation i x ( y ) refers to the x(y) component of the 
site i. 

Although the Hamiltonian we study is unconventional, it 
has all the relevant symmetries of the Hubbard model when 
J 3 = 1. In particular when /i = the Hamiltonian is invari- 
ant under the SU (2) spin and SU (2) pseudo-spin transforma- 
tions. When (j, ^ 0, the pseudo-spin symmetry is broken to the 
1/(1) fermion number symmetry. One can introduce repulsion 
or attraction by making U sufficiently positive or negative re- 
spectively. The important difference with the Hubbard model 
is that when U = the model is still strongly interacting and 
by setting J3 ^ 1 we can break the pseudo-spin symmetry. 
Further, the model simplifies in the U — > — 00 limit; in this 
limit the model can be mapped to the simple Hamiltonian 

H = ^J i -J, + (J 3 -l)JfJf. (6) 

Clearly, when J 3 = 1, one obtains the antiferromagnetic 
Heisenberg model involving pseudo-spins, while J 3 = leads 
to the XY model. When J3 = 1, U — ► 00 and /i^ = one gets 
the antiferromagnetic spin model 

H = J2Si-Si- ( 7 ) 

m 

An interesting aspect of this model is that the fermion sign 
problem can be solved using the meron cluster approach for 
< J3 < 1 when U < at any value of /Zj. Also when 
J 3 = 1 the sign problem can be solved when U > when /j, < 
U /2. Thus, we think the model offers a rich phase diagram 
and deserves to be investigated. In this paper we will consider 
— 00 < U < and investigate the physics when J3 = and 
J 3 = 1. We introduce disorder through //j = fi + where 
5 Hi is a random number distributed uniformly from — A to A. 

In order to probe superconductivity in this system we will 
focus on two observables. The first is the pair susceptibility 
defined as 

1 /T 1 JT 

*P = V £ / ^ l ^'^T' + P^r.jA (8) 

v Ui Jo Jo 
where V is the spatial volume, T is the temperature, and 

(9) 

is the pair correlation with P itJt , = P^p^^^y The 
second is the winding number susceptibility defined as 

1 jT 1 jT 

^w = ^E/ dT J o dr'[Cg jTt +C^ jTl ] (10) 
where 

Cg r ,(r) = ±Tr{e-^)»J,(i)e-^-T')B Mj)e -T>x} 

(11) 



is the current-current correlation function where J^(i) is the 
conserved fermion current, 



m = \ E [4 



(12) 



8=U 



at the site i in the fi = x,y direction. 

In order to estimate the importance of fermions we will also 
look at the density of singly occupied sites defined as 



— ^Tr[exp(-,&H"){n it +n tl -2n iT n a } 



and at the total density of electrons defined by 



E Tr I ex P(-^H n *T + ^1 } 
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(13) 



(14) 
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FIG. 2: An example of a bond configuration in one space and one 
time dimension. The configuration contains five clusters. 



where oj-p is the weight of the bond associated with the pla- 
quette V and takes one of the values given in Eq. J15i . 



A comparison of n s and n will tell us how many sites have 
formed local pairs. 



III. CLUSTER REPRESENTATION 

It is possible to rewrite the partition function of our model 
in discrete time in terms of a statistical mechanics of closed 
loops on a space-time latticed We first divide (3, the length 
in the Euclidean time direction, into M equal steps such that 
e = (3/M. Interactions between nearest neighbor sites are 
introduced in a checker-board type manner, so that on each 
time slice every site interacts with a unique neighbor. This 
then introduces 4 extra time slices for every e time step. In the 
cluster representation the nearest neighbor interactions occur 
in the form of three types of bond configurations on space- 
time plaquettes as shown in Fig. ^ Their weights are given 

by, 



L0 A 

u H = e' 
lu e = e' 



eJ 3 /4, e -eJ 3 /2 



- £ / 2 )/2 
e £ / 2 )/2 



(15) 



J 3 /4 (e -eJ 3 /2_ e -e/2 )/2 



Given a configuration {C} of bonds, one can connect them 
together to form many closed loops; we denote them by 
C a , a = 1, 2, N c . The partition function can then be writ- 
ten as 



Z = 



FIG. 1 : The three bond configurations on space-time plaquettes and 
their weights. The values of the weights are given in Eq. d 1 5 1 . 



fi(C a ) = 2cosh(-//cJ +a(C a )2ei 



(17) 



is the weight corresponding to each loop C a that arises due 
to the fermionic degrees of freedom associated with the loop. 
Here 



E 



(18) 



where the sum is over all all space-time points that belong to 
the cluster C a . Thus Sc a is just the size of the cluster. On the 
other hand 



E 

(iT)ec a 



(19) 



where w iT is +1 when the cluster is going forward in time 
and —1 when it is going backward in time at the site it. If 
the cluster moves horizontally, our convention is that the tem- 
poral direction is reversed. In order to determine Wi T one can 
start from any point and traverse the cluster in either direction. 
Note that 



W t {C a ) 



E 

(ir)eC Q 



(20) 



is the temporal winding of the cluster C a . Finally, the factor 
a(C a ) in Eq. d!7i is a sign factor associated with the cluster 
topology, that arises due to the fermion permutation signs. 22 
Following Ref. |22| we call the cluster a meron if a(C a ) — 
— 1. If Nh{C a ) is the number of horizontal hops in the cluster 
C a , then the cluster is a meron if and only if Nh(C a )/2 + 
Wt(C a ) is even. An example of a bond configuration in two 
dimensions is shown in Fig. [3] 

In the cluster representation it is easy to show that 



\P + ■ , + Pr . , 



2 cosh 



%r;jr')6C a ^ (21) 
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where 6<ir\jT')&C a imposes the restriction that both the space- 
time sites (ir) and (jr') belong to the cluster C a and 



{kT")GC a 



(22) 



where s l ^ J J is +1 while going from (ir) to (jr') and —1 
while continuing from (jr') to (ir). The winding number 
susceptibility is given by 



Xv 



41/ 

E 



W^ a cosh(|Mc a ) 
2fi(C Q ) 



W^gW^g* sinh(|^ Ca ) sinh(|^ Co) ) 



where Wu a refers to the spatial winding of the loop clusters 
along the spatial direction /i. The density of single occupation 
turns out to be 



— (E 



a{Cg)S Ca 2e^ s ^ 

4MF\^ n(c a ) 

and is a measure of the number of unpaired fermions. 
IV. DIRECTED-LOOP ALGORITHM 



A simple Monte Carlo algorithm for the current problem 
involves visiting every interaction plaquette and updating the 
bond configuration on it by replacing it with one of the three 
choices shown in Fig.^ Since the Boltzmann weight also de- 
pends on the structure of the loops formed by these bonds, the 
decision involves figuring out the connectivity of the sites of 
the plaquette (referred to as P,Q,R,S) due to the bonds on 
other plaquettes. This connectivity can be one of three types 
as shown in Fig. [5] Thus, choosing a new interaction involves 
finding weights of nine configurations (three bond configura- 
tions for three types of connectivity) and the new bond con- 
figuration can be found by using a heat bath or a Metropolis 
step. 

We have found that in the presence of a chemical poten- 
tial /x, this simple algorithm is inefficient. This behavior can 
be understood by noting that a chemical potential is similar 
to a magnetic field in a a quantum spin model and in the 
context of quantum spin models there is evidence that this 
type of naive algorithm becomes inefficient in the presence 
of magnetic fieldsi^&S Today it is well known that quantum 
spin systems in the presence of a magnetic field can be solved 



efficiently using the directed-loop algorithm. 25 However, un- 
til now this algorithm has been formulated only in the spin 
representation and not in the cluster representation. Unfortu- 
nately, the sign problem in the fermionic model can only be 
solved in the cluster representation. In this article we show 
how one can extend the directed-loop algorithm to the clus- 
ter representation which then leads to an efficient algorithm 
for the fermionic model even in the presence of a chemical 
potential. 

The basic idea behind the directed loop algorithm is to ex- 
tend the configuration space so that configurations that con- 
tribute to certain two point correlation functions (denoted 

( 2) 

{C iT j T i}) are sampled along with the configurations that con- 
\ tribute to the partition function (denoted {C}). The configu- 

/ rations Cj£j T , have two reference space-time points, ir and 

(23) ■? r '' during the directed-loop update one of these points, say 
ir, is held fixed while the other point jr' is moved around. 
The directed loop update begins with a configuration in the 
set {C} and chooses a site ir at random and probabilistically 

(2) 

creates a configuration in the set {C iT j T ,}, with ir = jr'. 
The probability of creating this configuration must satisfy de- 
tailed balance in order to produce configurations {C} and 

(2) 

{C^j. - T ;} with the correct Boltzmann weight. Once a config- 

(2) 

uration in the set {C iT ^ T ,} is created, the point jr' is moved 
around while satisfying detailed balance and thus sampling 

(2) 

other configurations in the set C iT - T , with the correct Boltz- 
mann weight. Finally, when the two points meet again, i.e.., 
when jr' = ir, the two points may be removed to obtain a 
configuration in the set {C} in accord with detailed balance. 
Thus, since at every step detailed balance is satisfied, it is easy 
to show that the directed-loop update, which starts from a con- 
figuration in {C} and ends on another configuration in {C}, 
satisfies detailed balance. During the loop update, all the sites 
encountered contribute to the two point correlation function. 

In the current work we have used the pair correlations to 
develop the directed loop algorithm. Thus, the configurations 
in {C\ T iT ,} are the ones that contribute to the pair correlations 
in Eq. 1211 . Thus, the weight of such a configuration is taken 
to be 



(24) 



FIG. 3: The three possible types of connectivity of the sites of a 
plaquette (P, Q, R, S) due to bonds on other plaquettes. 



2cosh(|4:f) 



{ II n(C a )Y[cjT,} (25) 



where the sites ir and jr' are forced to remain on the same 
cluster oq. The weight of a configuration in the set {C} is 



{n^ c «)n^}- 



(26) 



The essential steps of the directed loop algorithm are as fol- 
lows: 

(i) Start with the initial configuration which belongs to the 
set {C}. 

(ii) Select a space-time site ir at random and propose to 
create a configuration C\ T . T , assuming jr' = ir. Ac- 
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FIG. 4: Ten possible configurations for one of the three types of 
connectivity shown in Fig. [5] During the directed loop update, de- 
pending on the connectivity of the sites of the plaquette, a heat bath 
is used to pick one of ten possible configurations in order to move 
the site (jV). 



cept the proposal with probability 

( 2 cosh 



Min-j 



/ e ..('iT'jT') 



V(C ao 



(27) 



where C a<) is the cluster which contains the site it. 

(iii) If the proposal is not accepted then the update is com- 
plete. Otherwise we go on. 

(iv) The site jr' is moved to the next site by picking the 
plaquette it is connected to which is not the one just 
visited. Since each site is connected to two plaquettes 
the plaquette is unique unless one is at the beginning. 
In that case one chooses one of the two randomly. 

(v) Each plaquette update involves choosing one of ten 
possible configurations. These configurations depend 
on the connectivity of the plaquette and an example is 
given in Fig.0] One of these ten is chosen using a heat 
bath. 

(vi) Steps (iv)-(v) are repeated until jr' reaches it, at which 
stage the transition to the {C} sector is made with prob- 
ability 



2cosh(fA*£f /) 



,1 



(28) 



If the transition is made then the directed-loop algo- 
rithm ends, otherwise one goes back to step (iv) assum- 
ing one is in at the beginning of the loop. 

In the above algorithm, the pair susceptibility can be com- 
puted using the formula 
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(S) 



where S is the number of sites visited during the directed loop 
update. Other observables such as Xw and n s can be computed 
using the formula of Eq. i23\ and J24l> . 

We have tested the efficiency of the directed loop algorithm 
by comparing the results with spin model results in the limit 
U = — oo which can be obtained using the usual directed loop 
algorithm. Since each update of the plaquette in the fermion 
algorithm requires knowing the connectivity of the loop, the 
algorithm is indeed much slower than the directed loop algo- 
rithm of a spin model which does not require this step. Un- 
fortunately, this is a price one has to pay for being able to 
compute quantities in a fermionic theory. In Ref. 23 a trick 
was used to reduce the time to determine the connectivity of 
the clusters. The trick was to use a "tree" structure to store the 
information about the cluster connectivity which allowed one 
to obtain the relevant information in a time that grew like the 
logarithm of the cluster size. This was not implemented in the 
current work but could be implemented if necessary. 



V. RESULTS 

In this Section we discuss the results obtained from exten- 
sive simulations for lattice sizes up to L — 32. In our work 
we have fixed e = 0.25 in order to avoid changes in the time 
discretization errors. We have found that this value of e is rea- 
sonably small and the results at smaller values appear to join 
within our error bars. Further, since our desire is to understand 
universal physics of disorder, we believe that fixing e should 
not be a major concern since it only changes the transfer ma- 
trix by a small amount. For a given disorder realization we 
typically discard the first 1000 directed loop updates for equi- 
libration and then average over 20000 directed loop updates 
in blocks of 1000 to generate each of our statistics. All quan- 
tities plotted at a given value of A have been averaged over 20 
disorder realizations. 



A. Superconductivity with Disorder 

It is known from earlier studies 23 that our model has a low 
temperature superconducting phase when U = 0, J3 = and 
fi = A = 0. The superconductivity disappears at a finite tem- 
perature and the phase transition belongs to the expected BKT 
universality class. Here we study the effects of disorder on this 
system by keeping fj, = but A ^ 0. In the clean model the 
BKT predictions for the leading finite size scaling form of the 
pair susceptibility and the winding number susceptibility are 
known: 



Xp(L) 



AL 2- V t <T C 
A T>T C 



and 



Xw(L) 



B 



1 j i 

1 ^ 21og{L/L } 

Bexp(-L/L ) 



T~T C 
T>T C 



(30) 



(31) 



(29) 



where A, B and L are constants which depend on the tem- 
perature. We further expect < 77 < 0.25, and 77 = 0.25 with 
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FIG. 5: Pair susceptibility divided by L 1 ' 75 (xp/L 1 ' 75 , top, loga- 
rithmic scales) and winding number susceptibility (xw, bottom, lin- 
ear scales) as a function of lattice size for different magnitudes of 
disorder A at J 3 = 0, U = 0, fi = 0, and T = 0.25. The unfilled 
diamonds give the value of Xw at A = 1.65 in the XY limit obtained 
when U = — oo. The dashed horizontal lines are inserted as guides to 
the eye. The BKT transition appears to be near A ~ 1.5. 



B = 2 at the phase transitioni 2 ^ 2 ^ Here, instead of varying the 
temperature we fix the temperature at T = 0.25 and study the 
effects of disorder by increasing its strength through the pa- 
rameter A. If the BKT universality holds we then expect the 
same finite size behavior to be true where the constants A, B 
and Lq now depend on A. 

In Fig. [5] we plot our results for Xp and Xw as a function 
of L. The fits are shown in Table [I] The data appear to be 
consistent with a BKT transition around a critical disorder of 
~ 1.5. Note in particular the excellent power law behavior 
of Xp for A < 1.55. The expected form for Xw is seen for 
A < 1.35, but it does not fit very well for Aw 1.5. From the 
rapid decay of Xw (confirmed by saturation in x P , not shown), 



TABLE I: Fitting results for the pair susceptibility (x P ) and winding 
number susceptibility (Xw) for J3 = 0, U = 0, and /i = (Fig. [5}- 
The fitting formulas are Eqs. 1301 - 13 II : the \ 2 given is per degree of 
freedom from the fit. 
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0.198(9) 0.306 


2.49(7) 0.1(2) 0.409 


1.35 


0.46(1) 


0.198(7) 0.020 


2.14(3) 0.6(2) 0.251 


1.45 


0.474(8) 0.219(6) 0.834 


2.04(4) 0.7(3) 1.644 


1.55 


0.50(1) 


0.25(1) 0.676 


1.91(4) 0.8(4) 3.838 
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1.72(4) 1.4(4) 2.291 



the system is clearly no longer superconducting when A = 3. 

To obtain a good estimate of the critical disorder, A c , we 
assume that the forms fl30t-d3 li hold close to the transition 
and that the deviation of the exponent 77 and constant B is 
linear in A: 



r) = a(A - A c ) + 0.25 
B = 6(A-A c ) + 2.0. 



(32) 
(33) 



A joint fit of 77 and B to this form for A in the interval 
[1.35, 1.65] yields A c = 1.53(4). 

When U= — 00 our model reduces to the XY model in the 
pseudo-spin variables. As seen in Fig.|5J A = 1.65 is still in 
the superconducting phase in the XY limit. 29 Thus, the effect 
of fermions is to disorder the superconductor more quickly, as 
can be intuitively expected because of the increased entropy. 

As discussed earlier, when J3 = 1 and ^ = the model 
has an additional SU (2) pseudo-spin symmetry, as in the at- 
tractive Hubbard model. Thus, due to the Mermin-Wagner 
theorem, superconductivity is only possible when /i ^ 0. The 
J3 = 1 model also has been studied earlier in the absence of 
disorder, 24 and a BKT transition was established using uni- 
versal finite size scaling. We have extended these calcula- 
tions to the disordered regime. We again fix the temperature 
at T = 0.25 and study the effects of disorder with fx = \. 

Fig. Ogives our results and Table UTI shows the correspond- 
ing fits. The expected BKT transition behavior is indeed seen 
for A < 0.5. Note that the values obtained for 7/ are constant, 
within statistical error, in this range. Thus to extract the criti- 
cal disorder A c we use only the Xw data. A fit of the A < 0.5 
values to the form in Eq. J33l > yields A c = 0.13(5); the data 
for 77 is consistent with this value. The critical disorder found 
here is much smaller than in the J3 = case above, indicating 
as expected that superconductivity is weaker when J3 = 1. 



B. Fermionic degrees of freedom 

Since we are studying a strongly correlated electronic sys- 
tem with on-site attraction between spin-up and spin-down 
electrons, one might worry that the electrons have formed lo- 
cal pairs on each site and the model is essentially bosonic. 
Hence, it is important to look at observables that extract 
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FIG. 6: Pair susceptibility divided by L 1 ' 75 (Xp/L 1 - 75 , top, logarith- 
mic scales) and winding number susceptibility (Xw> bottom, linear 
scales) as a function of the lattice size for different magnitudes of 
disorder A at J 3 = 1, U = -0.1, fx = 1, and T = 0.25. The dashed 
horizontal lines are guides to the eye. Note the excellent BKT tran- 
sition behavior. 



fermionic information in the model. One such observable is 
the density of singly occupied sites defined in Eq. dl 31 . In 
Fig. we plot n s as a function of temperature when J3 = 
and /1 = 0. Since there is particle-hole symmetry at fi = we 
have 71 = 1. For this calculation we have fixed e = (3 /M = 0.25 
for (3 > 4 and M = 64 for (3 < 4. We show two different values 
of A - one for weak disorder and the other for strong disorder 
- at different values of U. The data shown was obtained at 
L = 16; we have observed that the density does not vary much 
as the lattice size increases. 

First we note that the density of singly occupied sites de- 
creases significantly as the temperature is lowered for U = 0. 
When the disorder is small the pairs begin to break close to the 



TABLE II: Fitting results for the pair susceptibility (x P ) an d me 
winding number susceptibility (xw) for J3 = 1, U= —0.1, and fi= 1 
(Fig.[6). The fitting formulas are Eqs. 1301 - 13 II ; the x 2 given is per 
degree of freedom from the fit. 





Xp 


Xw 


A 




B Lo x 2 


0.00 


0.365(8) 0.242(9) 0.081 


2.04(4) 0.5(2) 0.120 


0.25 


0.354(8) 0.235(9) 0.383 


1.97(3) 0.8(2) 0.277 


0.50 


0.346(7) 0.245(7) 0.460 


1.83(2) 1.0(2) 1.213 


0.75 


0.36(1) 0.29(1) 0.446 


1.58(2) 1.8(3) 1.952 


1.00 


0.35(1) 0.31(1) 0.083 


1.33(2) 2.5(3) 1.618 



critical temperature, while at large disorders the pairs break 
only at a temperature much higher than the critical tempera- 
ture. Further, as \U\ is reduced from infinity, the background 
density of singly occupied sites increases as expected; the in- 
crease is significant for T > 0.4, while it is modest for smaller 
T. The background density approaches a constant as T de- 
creases, leading us to conclude that fermionic excitations do 
exist for all values of T. On the other hand the change in 
disorder A has almost negligible effect on n s , especially near 
the phase transition (dotted lines in the figure). This leads 
us to conclude that the properties of the phase transition are 
most likely governed by a bosonic model such as the disor- 
dered quantum XY model, especially in the strongly disor- 
dered case. Thus, our work gives credence to the "dirty bo- 
son" scenario. 



co 0.2- 




(3(=1/T) 



FIG. 7: Density of singly occupied sites versus inverse temperature 
for A = 0.5 and A = 2 at J3 = 0, /1 = 0, and three different val- 
ues of U. The dotted line indicates the approximate critical value of 
j3 — 1 /T where the system undergoes a BKT transition to a super- 
conductor. The total particle density is n = 1 due to particle-hole 
symmetry. 
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VI. DIRECTIONS FOR THE FUTURE 

In this article we have studied the effects of disorder on a 
strongly correlated electronic system. Our model was known 
to be superconducting in the clean limit and our motivation in 
this work was to study the effects of disorder through a po- 
sition dependent chemical potential. Unfortunately, a naive 
extension of the meron cluster algorithm was found to be in- 
efficient in the presence of a disordered chemical potential. 
Hence in this work we constructed a new and efficient al- 
gorithm by combining the the meron-cluster approach with 
the directed-loop algorithm. Earlier work on the directed- 
loop algorithm involved quantum spin systems and always 
was constructed in the spin representation. In this work we 
have shown that the algorithm can be constructed even in the 
cluster representation, which is essential in the fermionic sys- 
tem in order to solve the fermion sign problem. Our new al- 
gorithm was quite successful and allowed us to compute the 
pair and winding number susceptibilities accurately. 

We found that disorder significantly suppresses supercon- 
ductivity and the system undergoes a phase transition which 
appears consistent with the BKT universality class. Although 
this scenario has been expected, 15 our work is the first, as far 
as we know, to study the universal scaling predictions of the 
BKT transition in a fermionic system with disorder. We could 
go to lattices as large as L = 32 thanks to our new algorithm. 



We found that when J3 — O superconductivity is stronger than 
when J3 = 1. 

We also found that there indeed are fermionic excitations 
in the system, but they are not affected by the disorder. The 
role of these fermions remains an interesting open question. 
For example do the background fermions form a Fermi-liquid 
in the weak disorder regime? If this is the case, then it would 
be interesting to ask whether the fermions become localized 
or do they remain extended. Is the phase transition between 
a superconductor and an insulator or whether it is a transition 
between a superconductor and a metal. Finally, although we 
have focused on attractive interactions in this work, we can 
study the repulsive model by setting U positive. In that case it 
is possible to add a chemical potential such that fj, < U/2 
when J3 = 1 without introducing a sign problem. 22 These 
studies have the potential to increase the fermionic effects fur- 
ther. 
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